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ABSTRACT 

An unstructured finite volume procedure has been developed for steady and transient 
thermo-fluid dynamic analysis of fluid systems and components. The procedure is 
applicable for a flow network consisting of pipes and various fittings where flow is 
assumed to be one dimensional. It can also be used to simulate flow in a component by 
modeling a multi-dimensional flow using the same numerical scheme. The flow domain 
is discretized into a number of interconnected control volumes located arbitrarily in 
space. The conservation equations for each control volume account for the transport of 
mass, momentum and entropy from the neighboring control volumes. In addition, they 
also include the sources of each conserved variable and time dependant terms. The source 
term of entropy equation contains entropy generation due to heat transfer and fluid 
friction. Thermodynamic properties are computed from the equation of state of a real 
fluid. The system of equations is solved by a hybrid numerical method which is a 
combination of simultaneous Newton-Raphson and successive substitution schemes. The 
paper also describes the application and verification of the procedure by comparing its 
predictions with the analytical and numerical solution of several benchmark problems. 

NOMENCLATURE 


A Area (in 2 ) 

Ao Pump Characteristic Curve Coefficient 

B 0 Pump Characteristic Curve Coefficient 

D Diameter (in) 

f Darcy Friction Factor 

g Gravitational Acceleration (ft/ sec 2 ) 

g c Conversion Constant (= 32.174 lb-ft/lb r sec 2 ) 

J Mechanical Equivalent of Heat (778 ft-lb/Btu) 

Kf Flow Resistance Coefficient (/lb r sec 2 /(lb-ft) 2 ) 

1^, Non-dimensional Rotating Flow Resistance Coefficient 
L Length (in) 

M Molecular Weight 

m Resident Mass (lb) 

• 

m Mass Flow Rate (lb/sec) 

p Pressure (lb/ in 2 ) 

Q Heat Source (Btu/sec) 



Re Reynolds Number (Re = puD/p.) 

R Gas Constant (lb r ft/lb-R) 

r Radius (in) 

S Momentum Source (lb f ) 

s Entropy (Btu/lbm -R) 

T Temperature (° F) 

u Velocity (ft/sec) 

V Volume (in 3 ) 

x Coordinate Direction 

y Coordinate Direction 

z Compressibility Factor 

Greek 

p Density (lb/ft 3 ) 

0 Angle Between Branch Flow Velocity Vector and Gravity Vector (deg). 
Angle Between Neighboring Branches for Computing Shear (deg) 
o> Angular Velocity (rad/sec) 

e Absolute Roughness (in) 

e/D Relative Roughness 

p Viscosity ( lb/ft-sec) 

y Specific Heat Ratio 

At Time Step (sec) 

t Time (sec) 


INTRODUCTION 

A numerical procedure capable of modeling fluid flow for both systems and components 
is of significant importance in aerospace and many other engineering industries. 
Although separate numerical methods are available for performing system and 
component analysis, there are not many attempts to develop an unified procedure 
applicable to both systems and components. Hardy Cross' developed a numerical method 
for calculating flow and pressure distribution in a flow network. This method has been 
widely used in modeling water distribution systems in municipalities. The Hardy Cross 
method assumes an equilibrium between pressure and friction forces in steady and 
incompressible flow. As a result it has not been successfully used in unsteady and 
compressible flow calculations where inertia force is important. 

Patankar and Spalding 2 developed a finite volume procedure to solve for Navier-Stokes 
equations in a structured co-ordinate system. Since the publication of the original paper in 
1972, there have been several developments reported to improve the numerical 
performance of the original algorithm. Several general purpose computer programs that 
have been developed based on these procedures, are increasingly being used for 
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component analysis. However, very few applications of the finite volume method have 
been reported for system analysis. Datta and Majumdar 3 used this solution algorithm to 
develop a calculation procedure for manifold flow systems. However, this procedure can 
not be used for simulating a flow network due to its dependence on a structured co- 
ordinate system. The purpose of this paper is to present a finite volume procedure based 
on an unstructured co-ordinate system that can be applied to analyze both systems and 
components. 

MATHEMATICAL FORMULATION 

All finite volume procedures require solution of mass, momentum and energy 
conservation equations in a flow domain consisting of interconnected control volumes. 
For turbulent flows, additional conservation equations of turbulence parameters are 
solved to calculate enhanced momentum exchanges due to turbulence. For reactive 
flows, specie conservation equations are solved in conjunction with other conservation 
equations. The main characteristic features of the present method are the following: 

Unstructured Co-ordinate System 

A flow domain is resolved into a network consisting of nodes and branches. The nodes 
are connected by branches. The nodes are classified into two categories: boundary and 
internal nodes. At a boundary node, pressure and temperature are known; at an internal 
node they are computed by solving conservation equations. Each internal node can be 
connected with an arbitrary number of internal and boundary nodes. At each node mass 
and entropy conservation equations are solved in conjunction with the equation of state. 
The momentum conservation equations are solved at all branches. This process of 
discretization allows the development of the set of conservation equations in an 
unstructured co-ordinate system. A schematic of a typical flow network is shown in 
Figure 1. 
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Figure 1 - Schematic of a Flow System Consisting of Nodes and Branches 
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Viscous Effect 

In conventional finite volume procedure, viscous effects are modeled by shear stress, 
which is expressed as a product of effective viscosity and local velocity gradient. In 
order to calculate the velocity gradient accurately, the finite volume procedure requires a 
large number of control volumes. It is not practical to use such a large mesh for a system 
level calculation. Therefore, in this method the frictional effect is modeled in the 

momentum equation by the following expression: 

. 2 

^ friction =K f m (1) 

where Kj is a function of Reynolds number, density, area and geometry of the branch in 
consideration. However for multi-dimensional flows, viscous effects are modeled by 
conventional shear stress method. 


Thermodynamic Formulation 

Entropy is calculated at every node using a second law analysis 4 of the control volume. 
The temperature and density at an internal node are calculated from pressure and entropy 
using a modified virial equation of state 5 . During each iteration cycle, thermodynamic 
properties are calculated to ensure thermodynamic equilibrium prevails during the course 
of attaining the solution. 

Solution Scheme 

The system of equations is solved by a hybrid numerical scheme which is a combination 
of the Newton-Raphson and successive substitution schemes. The equations which are 
strongly coupled are solved by the Newton-Raphson method. The equations which are 
weakly coupled are solved by the successive substitution scheme. The successive 
substitution scheme is also used to develop an initial guess for the Newton-Raphson 
scheme to ensure numerical stability. 

GOVERNING EQUATIONS 


Figure 2 displays a schematic showing adjacent nodes, their connecting branches, and the 
indexing system used in this paper. In order to solve for the unknown variables, mass and 
entropy conservation equations and the equation of the rmodynamic state are written for 
each internal node and flow rate equations are written for each branch. 


Mass Conservation Equation 


m t +AT ~ m. 

Ax 


j = n 
S mij 
j = l 


( 2 ) 
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Figure 2 - Schematic of Nodes, Branches and Indexing Practice 

Equation 2 requires that, for the transient formulation, the net mass flow from a given 
node must equate to the rate of change of mass in the control volume. In the steady state 
formulation, the left hand side of the equation is zero. This implies that the total mass 
flow rate into a node is equal to the total mass flow rate out of the node. 


Momentum Conservation Equation 

The flow rate in a branch is calculated from the momentum conservation equation 
(Equation 3) which represents the balance of fluid forces acting on a given branch. A 
typical branch configuration is shown in Figure 3. Inertia, pressure, gravity, friction and 
centrifugal forces are considered in the conservation equation. In addition to these five 
forces, a source term S has been provided to represent any external momentum sources. 
For example, this external momentum source term can be used to model a pump in a flow 
system. 
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Figure 3 - Schematic of a Branch Showing Gravity and Rotation 
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The two terms in the left hand side of the momentum equation represent the inertia of the 
fluid. The first one is the time dependent term and must be considered for unsteady 
calculations. The second term is significant when there is a large change in area or 
density from branch to branch. The first term in the right hand side of the momentum 
equation represents the pressure gradient in the branch. The pressures are located at the 
upstream and downstream face of a branch. The second term represents the effect of 
gravity. The gravity vector makes an angle (0) with the assumed flow direction vector. 
The third term represents the frictional effect. Friction was modeled as a product of Kf 
and the square of the flow rate and area. Kf is a function of the fluid density in the branch 
and the nature of the flow passage being modeled by tire branch. The calculation of Kf 
for different types of flow passages is described later in this paper. The fourth term in the 
momentum equation represents the effect of the centrifugal force. This term will be 
present only when the branch is rotating as shown in Figure 3. K TOt is the factor 
representing the fluid rotation. is unity when the fluid and the surrounding solid 
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surface rotates with the same speed. This term also requires a knowledge of the distances 
between the upstream and downstream faces of the branch from the axis of rotation. 


Entropy Conservation Equation 

The entropy conservation equation for node i, shown in Figure 2, can be expressed 
mathematically as shown in Equation 4. 
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The entropy generation rate due to fluid friction in a branch is expressed as 


S ijjen 


f 


y 

m ,j APij, viscous Kf \ 




PJX Pu T u J 


(4a) 


Equation 4 shows that for unsteady flow, the rate of increase of entropy in the control 
volume is equal to the rate of entropy transport into the control volume plus the rate of 
entropy generation in all upstream branches due to fluid friction plus the rate of entropy 
added to the control volume due to heat transfer. The MAX operator used in Equation 4 
is known as an upwind differencing scheme 2 which has been extensively employed in the 
numerical solution of Navier-Stokes equations in convection heat transfer and fluid flow 
applications. When the flow direction is not known, this operator allows the transport of 
entropy only from its upstream neighbor. In other words, the upstream neighbor 
influences its downstream neighbor but not vice versa. 

Equation of State 

A modified virial equation of state 5 is used to calculate the density from the computed 
pressure and temperature at a given node. 


i=N 


j=M 


P = z AiW + Z Bj(T)p 2j+l e- cp2,T > 


i=\ 




(5) 


A; (T) and Bj (T) are polynomials in T and 1/T 

T hi s equation was originally developed by Benedict, Webb and Rubin (BWR) and later 
modified by Bender 5 . This equation was the basis of the computer code GASP developed 
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by Hendricks et al 6 . This equation enabled PVT calculations to be made in the liquid and 
vapor phases. The derived properties of enthalpy and entropy were also obtained. The 
present method uses the GASP computer code to compute density and other 
thermophysical properties during iterative cycles of computation. 

Compressibility Factor 

For unsteady flow, resident mass in a control volume is calculated from the 
compressibility factor assuming a thermodynamic equilibrium. 


m = 


pV_ 

RTz 


( 6 ) 


The following table shows how each equation is used to calculate the unknown variables 
to demonstrate the problem of mathematical closure. 

Table 1 - Mathematical Closure 


Unknown 

variables 

Available equations to solve 

Pressure 

Mass Conservation Equation (Eqn. 1) 

Flow rate 

Momentum Conservation Equation (Eqn. 2) 

Temperature 

Second Law of Thermodynamics (Eqn. 3) 

Density 

Equation of State (Eqn. 4) 

Mass (Unsteady 
Flow only) 

Compressibility Factor (Eqn 5) 


Multi-Dimensional Flow 


Multi-dimensional conservation equations must account for the transport of mass, 
momentum and entropy into and out of the control volume from all directions in space. 
The scalar transport equations (i. e., mass and entropy conservation equations) can 
account for such transport because each internal node can be connected with multiple 
neighboring nodes located in space in any arbitrary location (Figure 2). Equations 2 and 
4 account for the transport of mass and entropy respectively from all directions in space. 
On the other hand, the momentum conservation equation (Equation 3) is one 
dimensional. Multi-dimensional momentum transport can be accounted for by 
incorporating two additional terms in the momentum equation. These terms include: a) 
momentum transport due to shear, and b) momentum transport due to the transverse 
component of velocity. 

These two terms can be identified in the two-dimensional steady state Navier-Stokes 
equation which can be expressed as: 
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The first term on the left hand side of Equation 7 corresponds to the inertia term in 
Equation 3. The second term on the left hand side of Equation 7 corresponds to the 
transverse momentum exchange. The first term on the right hand side of Equation 7 
corresponds to the pressure term. The second term on the right hand side of Equation 7 
corresponds to the gravity term. The third term on the right hand side of the equation is 
negligible (based on an order of magnitude argument). The fourth term on the right hand 
side corresponds to the friction term in the one dimensional momentum conservation 
equation. For multi-dimensional flow friction is modeled by shear force as discussed 
below. 


Momentum Transport Due to Shear 
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Figure 4 - Branch and Node Schematic for Shear Exchange 


Figure 4 represents a set of nodes and branches for which shear forces are exchanged. 
Let branch 1 2 represent the branch for which the shear force is to be calculated. Branches 
N12 and S12 represent the parallel branches which will be used to calculate the shear 
force on branch 12. Let YS be the distance between branches 12 and S12, and let YN be 
the distance between branches 12 and N12. Let AS be the shearing area between 
branches 12 and SI 2, while AN is the shearing area between branches 12 and N 12. 


The shear force on a control volume can be expressed as 


ft 




(Ax)(Ay)(Az) = p — — AxAz = p— A 
v ' Ay Ay 


shear 


( 8 ) 


The finite volume formulation for shear for the i th branch can be written as: 
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( 9 ) 


where i is the current branch, npi is the number of parallel branches to branch i, and ns, is 
the number of parallel solid walls to branch i. 

Transverse Momentum Transport 

The transverse momentum component of Equation 7 can be expressed in terms of a force 
per unit volume. 

5pvu (pvA u) 

~dy~ V * ( Ay ) ( A x X A yX Az ) = (A x)(A z)(pvAu) = ih^Au (10) 

Figure 5 represents a set of nodes and branches for which transverse momentum 
exchange will take place. Let branch 12 represent the current branch which will receive 
transverse momentum from the surrounding branches. Branch S12 represents an adjacent 
parallel branch, while branches SI and S2 represent the adjacent normal branches. 



Figure 5 - Branch and Node Schematic for Transverse Momentum Exchange 

The finite volume formulation for transverse momentum for the i* branch can be written 
as 

dpvu 
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f 


1 


max - X j rriijk cosGijk ),0|| 
V lc=l l n ij 


- (uij cos0ij^rnax|o, j] j^- ihijk cos0,,k| 


(ID 


Detailed derivation of equation (9) and (11) appear in Gl 'SSP User’s manual 7 . 

SOLUTION PROCEDURE 

In the sample circuit shown in Figure 1 , pressures and temperatures are to be calculated 
for the 5 internal nodes; flow rates are to be calculated in the 10 branches. The total 

10 

Draft paper for 37 th AIAA Aerospace Sciences Meeting 







number of equations can be evaluated from the following relationship: Number of 
equations = Number of internal nodes * Number of scalar transport equations + Number 
of branches. Therefore, the total number of equations to be solved is 20 (= 5 X 2 +10). 

There is no explicit equation for pressure. The pressures are implicitly computed from 
the mass conservation equation (Equation 2). The flow rates are calculated from 
Equation 3. The inertia and friction terms are nonlinear in Equation 3. The pressures and 
mass flow rates appear in the flow rate equations. The flow rates also appear in the 
entropy equation. The governing equations to be solved are strongly coupled and 
nonlinear and therefore they must be solved by an iterative method. 

There are two types of numerical methods available to solve a set of non-linear coupled 
algebraic equations: (1) the successive substitution method and (2) the Newton-Raphson 
method. In the successive substitution method, each equation is expressed explicitly to 
calculate one variable. The previously calculated variable is then substituted into the 
other equations to calculate another variable. In one iterative cycle each equation is 
visited. The iterative cycle is continued until the difference in the values of the variables 
in successive iterations becomes negligible. The advantages of the successive 
substitution method are its simplicity to program and its low code overhead. The main 
limitation, however, is finding an optimum order for visiting each equation in the model. 
This visiting order, which is called the information flow diagram, is crucial for 
convergence. Under-relaxation (partial substitution) of variables is often required to 
obtain numerical stability. 

In the Newton-Raphson method, the simultaneous solution of a set of non-linear 
equations is achieved through an iterative guess and correction procedure. Instead of 
solving for the variables directly, correction equations are constructed for all of the 
variables. The intent of the correction equations is to eliminate the error in each equation. 
The correction equations are constructed in two steps: (1) the residual errors in all of the 
equations are estimated and (2) the partial derivatives of all of the equations, with respect 
to each variable, are calculated. The correction equations are then solved by the Gaussian 
elimination method. These corrections are then applied to each variable, which 
completes one iteration cycle. These iterative cycles of calculations are repeated until the 
residual error in all of the equations is reduced to a specified limit. The Newton-Raphson 
method does not require an information flow diagram. Therefore, it has improved 
convergence characteristics. The main limitation to the Newton-Raphson method is its 
requirement for a large amount of computer memory. 

The present formulation employs a novel numerical scheme, SASS (Simultaneous 
Adjustment with Successive Substitution) which is a combination of the successive 
substitution and Newton-Raphson methods. The mass and momentum conservation 
equations are solved by the Newton-Raphson method. The entropy conservation 
equations are solved by the successive substitution method. The underlying principle for 
making such a division was that the equations which are more strongly coupled are 
solved by the Newton-Raphson method. The equations which are not strongly coupled 
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with the other set of equations are solved by the successive substitution method. Thus, 
the computer memory requirement can be significantly reduced while maintaining 
superior numerical convergence characteristics. To improve the convergence and 
stability of the numerical scheme, the successive substitution method is used to provide 
an initial guess for pressure and flow rate. Equations 2 and 3 are rewritten such that 
pressures and flow rates can be estimated at each node and branch. 

In a typical unsteady calculation, the SAAS procedure consists of the following steps: 

1. At the beginning of a new time step, provide the initial solution of all dependant 
variables in the flow domain e. g. pressure, resident mass, density and entropy at all 
internal and boundary nodes, flow rates at all branches. 

2. Begin the outer iteration loop; the purpose of this loop is to calculate entropy and 
density at all internal nodes and flow resistance in the branches. 

3. Solve mass conservation equation (Equation 2) in internal nodes, momentum 
conservation equation in branches (Equation 3) and equation of resident mass in 
internal nodes (Equation 6) by Newton Raphson method (Appendix - A ). 

4. Solve entropy conservation equation by successive substitution method (Appendix - 
B) 

5. Calculate density and temperature from the equation of state for calculated pressure 
and entropy at each internal node. Viscosity is also computed from the 
thermophysical property correlation for calculated pressure and temperature. 

6. Calculate flow resistance parameter (K*) of each branch. K f is a function of density 
and viscosity. 

7. Calculate the maximum difference in values of entropy, density and flow resistance 
parameters between successive iterations. Steps 3 to 7 constitute one iteration cycle. 

8. Repeat steps 3 to 7 until the maximum difference is less than the specified 
convergence criterion. Steps 2 to 8 constitute all operations required for one time step. 

9. Repeat steps 1 to 8 until final time is reached. 

COMPUTER PROGRAM 

The numerical procedure described in the previous sections was incorporated into a 
general purpose computer program named Generalized fluid System Simulation Program 
(GFSSP ) 1 . The computer program has three major parts. The first part consists of the 
subroutines for the preprocessor. The preprocessor allo ws the user to interactively create 
the flow network model consisting of nodes and branches. All of the input specifications, 
including the boundary conditions, are specified through the preprocessor. The second 
major part of the program consists of the subroutines that provide the initial conditions 
and then develop and solve all of the conservation equations in the flow network. The 
third part of the program consists of the thermodynamic property programs, GASP and 
WASP that provide the necessary thermodynamic arid thermophysical property data 
required to solve the resulting system of equations. 
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Figure 6 shows GFSSP’s process flow diagram. The user runs the interactive 
preprocessor to generate the input data file. The input data file contains all the 
information necessary for the model. The solver module reads the input data file and 
produces the solution in conjunction with the thermodynamic property programs. 



Figure 6 - GFSSP Process Flow Diagram 
RESULTS 

The feasibility, robustness and verification of the proposed method have been 
demonstrated by simulating five example problems. They are: 

1 . Flow system consisting of a pump, valve and pipe line. 

2. Water distribution network. 

3. Compressible flow through a converging-diverging nozzle. 

4. Blow down of a pressurized tank. 

5. Recirculating flow in a square cavity 

Example 1 - Flow System Consisting of a Pump, Valve and a Pipe Line 

A problem commonly encountered in fluid engineering is to match a pump’s 
characteristics with the operating system’s characteristics. The designer needs to know 
the flow rate in the system and the power consumed by the pump. The system considered 
for this example is shown in Figure 7. 

The fluid system shown in Figure 7 was simulated with a finite volume model consisting 
of four nodes and three branches as shown in Figure 8. 
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Figure 7 - Flow system consisting of a pump, valve and pipe line. 
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Figure 8 - Finite volume Model of Flow System Described in Example 1. 


Pump Model 


The pump was modeled as a momentum source in branch 12. The pump characteristics 
were expressed as: 


Where: 


. 2 

Ap = A 0 + B 0 m (12) 


A p = Pressure rise, lbf/ft 2 
m = Flow rate, lbm/sec 

The momentum source, S, in Equation 3 was then expressed as: 
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S = Ap A 


Valve Model 


The resistance in branch 23 representing the valve was computed by the two-K method 
of Hooper ®. For this option, Kf of Equation 3 was expressed as: 


K f = 


Ki / Re+ Koo(l + 1 / D) 
^§cPu A 


(13) 


Where: 

Kj= K for the fitting at Re =1 ; 

K oo = K for the fitting at Re = oo ; 

D = Internal diameter of attached pipe, in. 

The constants K, and Koo for common fittings and valves are listed in reference 8. 

Pipe Model 

The resistance in branch 34 was computed from the friction factor in the pipe line. The 
resistance coefficient, Kf for a pipe with length, L, diameter, D, and surface roughness, e 
was expressed as: 


K f = 


8fL 


2 5 

Pu 71 D 8 C 


(14) 


Where p u is the density of the fluid at the upstream node of a given branch. 


The Darcy friction factor, f, is determined from the Colebrook Equation 9 which is 
expressed as: 


47 


= -2 log 


s 2.51 


3.7D Rejf 


(15) 


Where zfD and Re are the surface roughness factor and Reynolds number respectively. 


It took 28 iterations to satisfy the convergence criterion. The solutions are shown in 
Tables 2 and 3. Table 2 shows pressures, temperatures, compressibility factors and 
density in nodes 2 and 3. It may be noted that water was treated as a compressible fluid 
and the compressibility factor is very low as expected. 
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Table 2 . Predicted Solution of Example 1 at Internal Nodes 


INTERNAL 
NODE 

2 
3 

Table 3. Predicted Solution of Example 1 at Branches 


NODES 

P{PSI) 

0 . 22 90E+03 
0 . 22 88E+03 


T (F) 

0 .6003E+02 
0 .6003E+02 


Z 

0 . 1186E-01 
0 . 1185E-01 


RHO 

(LBM/FT"3) 
0 .6241E+02 
0 .6241E+02 


BRANCHES 


BRANCH KFACTOR 

DELP 

FLOW RATE 

VELOCITY 

REYN. NO. 

MACH NO. 

{LBF-S A 2 / (LBM-FT) A 2) 

(PSI) 

(LBM/SEC) 

(FT/SEC) 



12 0 . 000E+00 

-0 . 214E+03 

0 . 191E+03 

0 . 2 I9E+01 

0 . 241E+06 

0 . 183E-02 

23 0.764E-03 

0 . 193E+00 

0 . 191E+03 

0 . 1 >6E+02 

0 . 644E+06 

0 . 13 0E- 01 

34 0 . 591E+00 

0 . 214E+03 

0 . 191E+03 

0 . 156E+02 

0 . 644E+06 

0.130E-01 


Table 2 shows Kf , pressure drop, flow rate, velocity, Reynolds number and Mach number 
in each branch. The predicted flow rate is 191 lbm/sec and the pump is supplying 214 
psid pressure rise to meet the system requirement. The predicted system and pump 
characteristics are shown in Figure 8 which also provides verification of the predicted 
operating point shown in Table 3. 



0 500 1000 1500 2000 2500 3000 

Mas* Flow Rat* (Ibm/*) 

Figure 8 - Fluid System Operating Point 
Example 2 - Water Distribution Network 

In Example 1 we analyzed a single line pipe flow problem commonly encountered by 
pipeline designers. In this example, we consider an example associated with multipath 
systems which are commonly known as flow networks. In general, water supply systems 
are considered as flow networks, since nearly all such systems consist of many 
interconnecting pipes. A ten pipe (commercial steel) distribution system is shown in 
Figure 9. 
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Figure 9 - A Schematic of Water Distribution Network 
Table 4 - Water Distribution Network Branch Information 


Branch 

Length(inches) 

Diameter(inches) 

Roughness 

Factor 

12 

120 

6 

0.0018 

25 

2400 

6 

0.0018 

27 

2400 

5 

0.0018 

57 

1440 

4 

0.0018 

53 

120 

5 

0.0018 

56 

2400 

4 

0.0018 

64 

120 

4 

0.0018 

68 

1440 

4 

0.0018 

78 

2400 

4 

0.0018 

89 

120 

5 

0.0018 


The system shown in Figure 9 is modeled by using 9 nodes and 10 branches as shown in 
Figure 1. The fluid was assumed incompressible. Nodes 1, 3, 4 and 9 are boundary nodes 
where the pressures are prescribed. Node 1 represents the inlet boundary node. Nodes 3, 
4 and 9 are outlet boundary nodes. All of the remaining nodes (2, 5, 6, 7 and 8) are 
internal nodes where the pressures are calculated. All of the branches in this circuit 
simulate pipes. The length, diameter and roughness factor of all branches are given in 
Table 4. 

Figure 10 shows a comparison of flow rates between the present predictions and the 
Hardy Cross method. The comparison appears reasonable considering the fact that the 
Hardy Cross method assumes a constant friction factor in the branch while the present 
method computes the friction factor for each branch during every iteration. Therefore, as 
the flowrates change the friction factor also changes. 
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Figure 10 - A Flow Rate Comparison Between GFSSP and Hardy Cross Method 
Predictions 

Example 3 - Compressible Flow in a Converging-Diverging Nozzle 

In this example we will consider compressible flow in a converging-diverging nozzle to 
demonstrate the method’s capability to handle compressibility. One of the characteristics 
of compressible flow in a duct is that the flow rate becomes independent of exit pressure 
after reaching a threshold flow rate. This threshold value is known as the choked flow 
rate and it is a function of inlet pressure and temperature. Flow in a confined duct 
becomes choked when the flow velocity equals the local velocity of sound. The purpose 
of this example is to investigate how accurately the present procedure can predict the 
choked flow rate in a converging-diverging nozzle. 

The converging-diverging nozzle considered for this example is shown in Figure 1 1 . The 
fluid considered was steam at 150 psia and 1000 °F. The nozzle back pressure was varied 
from 134 psia to 15 psia. 


T~ 

0.492 in. 


0.158 in. - » 


Figure 1 1 - Converging-Diverging Steam Nozzle Schematic 

The fluid system shown in Figure 11 was simulated by seventeen nodes and sixteen 
branches. Nodes 1 and 17 are the boundary nodes representing the inlet and outlet of the 
nozzle. All of the remaining nodes are internal nodes connected in series. The outlet 
boundary node pressures were varied to include 134, 100, 60, 30 and 15 psia. Table 5 lists 
the model predicted mass flow rates with varying exit pressures. As expected, the mass 
flow rate increased as the exit pressure was decreased until the pressure ratio decreased 
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below the critical pressure ratio. At this point and below, the mass flow rate remained 
constant due to choking of the flow at the nozzle throat. 

Table 5 - Predicted Mass Flow Rate with Varying Exit Pressure 


p 

1 exit 

(psia) 

m (lbm/s) 

134 

0.258 

100 

0.301 

60 

0.308 

30 

0.308 

15 

0.308 


The isentropic flow rate at the choked point was calculated to be 0.303 lbm/sec from the 
following relation: 


m = A 


throat inlet 


CTQ 

o 

' 2 ] 

(K-3 

R Tintet > 

^ Y + 1/ 



Example 4 - Simulation of the Blow Down of a Pressurized Tank 


(16) 


In the previous examples we considered the simulation of steady state flow in a given 
flow circuit. In this example we will employ the capabilities of the unsteady flow 
formulation to simulate the process of blowing down a pressurized tank. 

Consider a tank with an internal volume of 10 containing nitrogen gas at a pressure 
and temperature of 100 psia and 80 °F respectively. The nitrogen is discharged into the 
atmosphere through an orifice with a 0. 1 inch diameter until the pressure in the tank drops 
to 50 psia. The purpose of this example is to determine the blow down time and the 
pressure, mass flow rate and temperature history of the isentropic blow down process. 
These predicted values will then be compared with the analytical solution. 


/ 1 


Tank 


v = to ft 3 


Initial Conditions 

p { = 100 psia 
T, = 80°F 


d = 0.1 in 


Atmosphere 
p = 14.7 psia 



Figure 12 - Physical Schematic (a) and GFSSP Model (b) for Venting Nitrogen Tank 

19 

Draft paper for 37 th AIAA Aerospace Sciences Meeting 




The physical schematic for Example 8 is shown in Figure 12(a) and a schematic of the 
corresponding finite volume model is shown in Figure 12(b). The venting process can be 
modeled with two nodes and one branch. Node 1 is an internal node which represents the 
tank. 


Analytical Solution : 

The differential equation governing an isentropic blow down process can be written as: 


N(l-3r)/2y 
- 
1 Pj 


diplp.) y A , -( 2 ) 

ST = 


(17) 


This is an initial value problem and the initial conditions are: 


T = 0, — = 1 


The analytical solution for p / Pj is given by Moody 1 as 


Pi 


1 + 


Y -1 


2 ) 


,(r+l)/2(y-l) 


vy +b 


Y ScPj At 
I P , V 


-2y /(y-l) 


(18) 


The analytical and finite volume solutions obtained by GFSSP are compared in Figure 
13. Excellent agreement was observed between the two solutions. 



Figure 13 - Comparison of the Predicted Pressure History by Finite Volume Method 
using GFSSP and the Analytical Solution 
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Example 5 - Recirculating Flow in a Square Cavity 


Flow in a square cavity is induced by moving its top wall at a constant speed as shown in 
Figure 14. The density of the fluid is assumed constant at 1.00 lb m /ft 3 , and the viscosity of 
the fluid is assumed to be 1.00 lb m /(ft sec). The bottom and side walls are fixed. The top 
wall is moving to the right at a constant speed of 1 00 ft/sec. The corresponding Reynolds 
number for this situation is Re = 100. 


Due to the complexity of the flow field, an analytical solution of this situation is not 
practical. Instead of an analytical solution, a well known numerical solution by Odus 
Burggraf" was used as the benchmark solution. Burggraf used a 51x51 grid in his model 
of the square cavity. 



Figure 14 - Flow in a Shear Driven Square Cavity 

The finite volume model of the driven cavity consists of 49 nodes (48 of which are 
internal) and 84 branches. For numerical stability, one node (Node 1) was assigned to be 
a boundary node with an arbitrary pressure of 100 psi. A unit depth (1 inch) was assumed 
for the required areas. The model is shown schematically in Figure 15. Modeling details 
are provided in GFSSP user’s manual 7 
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Figure 15 - Finite Volume Model of Flow in a Shear Driven Square Cavity 

Figure 16 shows a comparison between the benchmark numerical solution and the finite 
volume predictions along a vertical plane at the horizontal midpoint. As can be seen in 
Figure 16, the results of this coarse grid model compare very favorably with the 
benchmark numerical solution of Burggraf. 



Dimensionless Velocity 


Figure 16 - Shear Driven Square Cavity Centerline Velocity Distribution 
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CONCLUSIONS 


1 . The present paper has described a generally applicable numerical method to perform 
thermo-fluid dynamic analysis of a fluid system or component. 

2. A novel unstructured co-ordinate system has been used; this allows the development 
of mass, momentum and entropy conservation equations in a complex flow network 
defined by interconnected nodes and branches. 

3. The method described here uses a generalized momentum equation which considers 
an one-dimensional form for system level calculations. Fluid friction is calculated 
from the friction factor or loss coefficients. For component level analysis, the multi- 
dimensional form of the momentum equation is used. Fluid friction is, then, 
calculated from the local velocity gradient and viscosity. In addition, transverse 
transport of momentum is also computed. 

4. Thermodynamics and fluid dynamics of the flow is linked through entropy and 
pressure. The second law formulation allows the calculation of entropy generation 
through heat and fluid flow. All thermodynamic and thermophysical properties are 
computed from local values of pressure and entropy. 

5. The system of equations is solved by a novel numerical method which is a 
combination of successive substitution and Newton Raphson methods. The method 
has demonstrated a remarkable stability for both system and component level 
analysis. 

6. The method has been incorporated into a general purpose computer program which is 
extensively being used in many aerospace engineering applications and the computer 
program, GFSSP, is available through COSMIC, NASA’s Software Technology 
Transfer Center. 
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APPENDIX - A 

Newton-Raphson Method of Solving Coupled Nonlinear System of 

Algebraic Equations 


The application of the Newton-Raphson Method involves the following 7 steps: 
1 . Develop the governing equations. 

The equations are expressed in the following form: 

/l(*j >*2>*3> x n ) = 0 

■^2^ X \ >* 2 ’* 3 ’ x n ) = 0 

f n^ x i’*2’*3> xn) = 0 

If there are n number of unknown variables, there are n number of equations. 
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2. Guess a solution for the equations. 


Guess x\ ,X 2 ,^ 3 , x n asan initial solution for the governing equations 

3. Calculate the residuals of each equation. 

When the guessed solutions are substituted into Equation A-l, the right hand side of the 
equation is not zero. The non-zero value is the residual. 

f l( x p*2’*3’ *«) = R ] 

f 2 ( x \> x 2’ x 3> x *n) = R 2 (A _ 2 ^ 


fnl** V x 2,4. x n) = R n 

The intent of the solution scheme is to correct x\,x2,x^, x n with a set of 

corrections jq,jc2>*3> x n such that R\,R2, R 3, ,R n are zero. 

4. Develop a set of correction equations for all variables. 

First construct the matrix of influence coefficients: 

3/j df\ 9/j 5/j 

d x\ 3*2 & x 3 9 x n 

df 2 2 2 2 

dx\ d X2 dx 2 dx n 


d _£n d _£n d _£n d _£n 

3x] 9^2 3x3 9 x n 

Then construct the set of simultaneous equations for corrections: 


Sfl 

3/j 3/j 

3/i - 



3 x] 

9x2 ^ x 3 

d x n 

xj 


*1 

9/ 2 

9/2 5/2 

3/ 2 

*2 


*2 

3x] 

9x2 ^ x 3 

9 


- 



5 At 5 -A 





3xj 

3x2 ^ x 3 

3x„_ 
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5. Solve for jq ,%2 ,*' 3 , x' n by solving the simultaneous equations. 

6 . Apply correction to each variable. 

7. Iterate until the corrections become very small. 


APPENDIX - B 

Successive Substitution Method of Solving Coupled Nonlinear System of 

Algebraic Equations 

The application of the successive substitution method involves the following steps: 

1 . Develop the governing equations: 

*i = /l(*j,x 2>*3> x n) 

*2 = ^2 (j V*2’ x 3> *n) rR _ n 


X„ - f n {x^,x 2 ,^ 3 , x n ) 

If there are n number of unknown variables, there are n number of equations. 

2. Guess a solution for the equations: 

««« » 

Guess x \ , *2 > X2 , x n as an initial solution for the governing equations. 

3. Compute new values of x 1 .x 2 .x 3 , x n by subsi ituting x 1 .x 2 .x 3 , x* n in the 

right hand side of Equation B-l. 

4. Under-relax the computed new value: 
x = (l -a)jc* + ax 

where a is the under-relaxation parameter. 

jjc j|( )|{ 

5. Replace jci ,X 2 ,X 3 , x n with the computed value of x 1 .x 2 .x 3 , x n from 

Step 4. 

6 . Repeat Steps 3 to 5 until convergence. 
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